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Abstract 

We solve the geodesic deviation equations for the orbital motions in 
the Schwarzschild metric which are close to a circular orbit. It turns out 
that in this particular case the equations reduce to a linear system, which 
after diagonalization describes just a collection of harmonic oscillators, with 
two characteristic frequencies. The new geodesic obtained by adding this 
solution to the circular one, describes not only the linear approximation of 
Kepler's laws, but gives also the right value of the perihelion advance (in the 
limit of almost circular orbits) . We derive also the equations for higher-order 
deviations and show how these equations lead to better approximations, 
including the non-linear effects. The approximate orbital solutions are then 
inserted into the quadrupole formula to estimate the gravitational radiation 
from non-circular orbits. 
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1 Introduction 



The problem of motion of planets in General Relativity, considered as test par- 
ticles moving along geodesic lines in the metric of Schwarzschild's solution, has 
been solved in an approximate way by Einstein who found that the perihelion 
advance during one revolution is given in the near-Keplerian limit by the formula 

6ttGM 
a(l - e 2 



where G is Newton's gravitational constant, M the mass of the central body, a 
the greater half- axis of planet's orbit and e its eccentricity. 

This formula is deduced from the exact solution of the General Relativistic 
problem of motion of a test particle in the field of Schwarzschild metric, which 
leads to the expression of the angular variable <p as an elliptic integral, which is 
then evaluated after expansion of the integrand in terms of powers of the small 
quantity 

The formula has been successfully confronted with observation, giving excellent 
fits not only for the orbits with small eccentricities (e.g., one of the highest values 
of e displayed by the orbit of Mercury, is e = 0.2056 ), but also in the case when e 
is very high, as for the asteroid Icarus (e = 0.827), and represents one of the best 
confirmations of Einstein's theory of gravitation. In the case of small eccentricities 
the formula (JT|) can be developed into a power series: 

A(j) = (l + e 2 + e 4 + e 6 + ...). (2) 

a 

One can note at this point that even for the case of planet Mercury, the series 
truncated at the second term, i.e., taking into account only the factor (1 + e 2 ) will 
lead to the result that differs only by 0.18% from the result predicted by relation 
(|T|), which is below the actual error bar. 

This is why we think it is useful to present an alternative way of treating 
this problem, which is based on the use of geodesic deviation equations of first 
and higher orders. Instead of developing the exact formulae of motion in terms 
of powers of the parameter we propose to start with an exact solution of a 
particularly simple form (i.e. a circular orbit with uniform angular velocity), and 
then generate the approximate solutions as geodesies being close to this orbit. 

One of the advantages of this method is the fact that it amounts to treating 
consecutively systems of linear equations with constant coefficients, all of them 
being of harmonic oscillator type, eventually with an extra right-hand side being 
a known periodic function of the proper time. The approximate solution obtained 
in this manner has the form of a Fourier series and represents the closed orbit as a 
superposition of epicycles with diminishing amplitude as their circular frequencies 
grow as multiples of the basic one. This approach is particularly well-suited for 
using numerical computations. An example is provided by the computation of 
gravitational radiation from non-circular orbits, for which we use the well-known 
quadrupole formula. 
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2 Geodesic deviations of first and higher orders 

Of many equivalent derivations of geodesic deviation equation we present the one 
which most directly leads to the results used in subsequent applications. Given a 
(pseudo)-Riemannian manifold V4 with the line element defined by metric tensor 

ds 2 = (x x ) dx»dx v , (3) 

a smooth curve x x (s) parametrized with its own length parameter (or proper time) 
s is a geodesic if its tangent vector u p = (dx^/ds) satisfies the equation: 

u x v xU » = o & ^ = ^r + r V A ^ = - ( 4 ) 

where denote the Christoffel connection coefficients of the metric g^. 

Suppose that a smooth congruence of geodesies is given, of which the geodesies 
are labeled by a continuous parameter p: x M = x p (s,p), such that the two inde- 
pendent tangent vector fields are defined by: 

M ^( S) p) = _ and n^(s,p) = —. (5) 

It is easily established that the rates of change of the tangent vectors in the mu- 
tually defined directions are equal: 

Ay ^ A ^ M ^ Du^ Dn^ d 2 x fl ^ ^ M dx x dx p 

Dp Ds dpds xp dp ds 

by virtue of the symmetry of Christoffel symbols in their lower indices. 

The Riemann tensor can be defined using covariant derivations along the two 
independent directions of the congruence: 



[ M A V A ,r/V p ] 3" 



D D D D 

Ds Dp Dp Ds 



Y H »° ds dp {f) 



Replacing Y p by in the above formula, we get 

[u x V x , n'Vj u» = iV -^ ^ = iV u x u° n?. (8) 

By virtue of the geodesic equation (H) and Eq. (||), this can be written as 

n n 7/ M D 2 r> p 

u x V x (n'V/) = Ws^p-=^ = Rx ^ U U ° n ^ (9) 

This first-order geodesic deviation equation is often called the Jacobi equation, and 
is manifestly covariant. 
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In certain applications, Eq. (pD can be replaced by its more explicit, although 
non-manifestly covariant version: 



d 2 n p 



dn p 



_ + 2 r;y- + a „r; pU W 



0. 



(10) 



In this form of the geodesic deviation equation one easily identifies the relativistic 
generalizations of the Coriolis-type and centrifugal-type inertial forces, represented 
respectively by the second and third terms of Eq. (|T0| ) . 

The geodesic deviation can be used to construct geodesies x p (s) close to a 
given reference geodesic Xq(s), by an iterative procedure as follows. Let the two 
geodesies be members of a congruence as above, with 



x fJ '(s) = x p (s,p), a 
It follows by direct Taylor expansion, that 

dx p 



Xn(s) 



a^(s,Po)- 



x p (s,p) = x fl (s } po) + (p-p ) 



dp 



1 , , 2 <9V 



+ 



(*:P0) 



+ 5x"(s) + i <SV(s) + j 5 3 x»(s) + 



(12) 



where the more compact notation 8 n x tl (s) describes the n* h -order geodesic devia- 
tion. Because (p — po) is supposed to be a small quantity, for convenience we may 
denote it e. The first-order deviation is a vector, &c M (s) = (p — po) n oi s ) 
But the second-order deviation is not a vector, and is given by 



en r n 



S z x 



(p-p ) 2 (If - T>V) = e 2 (if - I*nV) 

where the covariant second-order deviation vector b p is defined by 

Drf drf 



Dp dp 



(13) 



(14) 



Straightforward covariant differentiation of Eq. @, plus use of the Bianchi and 
Ricci identities for the Riemann tensor, implies that this second-order deviation 
vector 6 M (s) satisfies an inhomogeneous extension of the first-order geodesic devi- 
ation equation: 



r> 2 h p 



Xpcr 



Dn u 



■ (15) 



A more detailed formal derivation of this equation is given in appendix 2. 

A rigorous mathematical study of geodesic deviations up to the second-order, 
as well as geometric interpretation, but using a different derivation, was presented 
in Ref. [0. Also, a Hamilton- Jacobi formalism has been derived in Refs. 0, which 
was applied to the problem of free falling particles in the Schwarzschild space-time 
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[|J . Fine effects resulting from the analysis of geodesic deviations of test particles 
suspended in hollow spherical satellites have been discussed in Ref. ||. 

Obviously the procedure can be extended to arbitrarily high order geodesic 
deviations 5 n x' M (s). This is of considerable practical importance, as it allows to 
construct a desired set of geodesies in the neighborhood of the reference Xq (s), when 
the congruence of geodesies is not given a priori in closed form. Indeed, all that 
is needed is the set of deviation vectors (tIq(s), 6q(s), ...) on the reference geodesic; 
obviously these vectors are completely specified as functions of s by solving the 
geodesic deviation equations (f|), (|i~5|) and their extensions to higher order, for 
given Xq(s). 

As in the case of the first-order deviation, it is sometimes convenient to write 
equation ([15]) in the equivalent but non-manifest covariant form 

d 2 h^ dh a rln u 

^ + d p i* « V&" + 2 ri ?L = 4 (d x r» p + k p ry (« V - ,V) 

+ (T w d T T% + 2r£ T d p Vl v - d v d a ry (t*VnV - uVnV). 

(16) 

An equation for the 3rd-order deviation is presented in Appendix 1. 

The non-manifestly covariant geodesic deviation equations are often better 
adapted to deriving successive approximations for geodesies close to the initial 
one. Starting from a given geodesic x^ (s) we can solve Eq. fllCf ) and find the first- 
order deviation vector (s). Then, inserting (s) and n M (s), by now completely 
determined, into the system (|16|), we can solve and find the second-order deviation 
vector 6 M (s), and subsequently for the true second-order coordinate deviation S 2 x^, 
and so forth. As an example, below we describe non-circular motion, along with 
Kepler's laws (in an approximate version), together with the relativistic perihelion 
advance, starting from a circular orbit in Schwarzschild metric. 

Although for orbital motion in a Schwarzschild background we have at our 
disposal the exact solutions in terms of quadratures (with integrals of elliptic or 
Jacobi type), our approach is particularly well-suited for numerical computations, 
because in appropriate (Gaussian) coordinates the geodesic curves can display a 
very simple parametric form, and all the components of the 4-velocity and other 
quantities reduce to constants when restricted to that geodesic. 

In this case equation fllCf ) reduces to a linear system with constant coefficients, 
which after diagonalization becomes a collection of harmonic oscillators, and all 
that remains is to find the characteristic frequencies. In the next step, we get a 
collection of harmonic oscillators excited by external periodic forces represented 
by the right-hand side of ([16]), which can also be solved very easily, and so forth. 

In the third order, the presence of resonances giving rise to secular terms could 
in principle lead to instability of the orbit we started with; but this phenomenon 
can be dealt with by Poincare's method ||, according to which such terms can be 
eliminated if we admit that the frequency of the resulting solution is also slightly 
modified by the exterior perturbation, and can be expanded in a formal series in 
successive powers of the initial (small) deformation parameter. 
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At the end, the deviation becomes a series of powers of a small parameter 
containing linear combination of characteristic frequencies appearing on the right- 
hand side, which are entire multiples of the basic frequency, also slightly deformed. 
This description of planetary motion as a superposition of different harmonic mo- 
tions has been first introduced by Ptolemaeus in the II century M. We shall now 
analyse the simplest case of circular orbits in Schwarzschild geometry. 



3 Circular orbits in Schwarzschild metric 

Let us consider the geodesic deviation equation starting with a circular orbit in 
the field of a spherically-symmetric massive body, i.e. in the Schwarzschild metric. 
The circular orbits and their stability have been analyzed and studied in several 
papers [§, [|, [10[ and books, e.g. the well-known monograph by Chandrasekhar 





The gravitational field is described by the line-element (in natural coordinates 
with c = 1 and G = 1) 



g^dx^dx" = -ds 2 = -B{r)dt 2 + dr 2 + r 2 (d6 2 + sin 2 9 dcj) 2 ) , (17) 

2M , , 



B(r, 

with 

B(r) = 1 - 

r 

We recall the essential features of the solution of the geodesic equations for a test 
particle of mass m « M. As the spherical symmetry guarantees conservation of 
angular momentum, the particle orbits are always confined to an equatorial plane, 
which we choose to be the plane 9 = tt/2. The angular momentum J is then 
directed along the z-axis. Denoting its magnitude per unit of mass by I = J/m, 
we have 

T = 4- I") 

ds r z 

In addition, as the metric is static outside the horizon r + = 2M, it allows a time- 
like Killing vector which guarantees the existence of a conserved world-line energy 
(per unit of mass m) e, such that 

dt 



ds 1 — 



2M 



(20) 



Finally, the equation for the radial coordinate r can be integrated owing to the 
conservation of the world-line Hamiltonian, i.e. the conservation of the absolute 
four-velocity: 

as J \ r J \ r A 

From this we derive a simplified expression for the radial acceleration: 

d 2 r M ( i 2 \ ( 3ikf 



1-— 1 + 3 ■ (21) 
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The equation (|2T| ) can in principle be integrated directly; indeed, the orbital func- 
tion r(<ft) is given by an elliptic integral fT2| , |13| . However, to get directly an 
approximate parametric solution to the equations of motion one can also study 
perturbations of special simple orbits. In the following we study the problem for 
bound orbits by considering the first and second-order geodesic deviation equa- 
tions for the special case of world lines close to circular orbits. 

Observe that for circular orbits r = R = constant, the expressions for dr/ds, 
Eq. (pTD, and d 2 r/ds 2 , Eq. fl22|), must both vanish at all times. This produces two 
relations between the three dynamical quantities (R,e,£), showing that the circu- 
lar orbits are characterized completely by specifying either the radial coordinate, 
or the energy, or the angular momentum of the planet. In particular, the equation 
for null radial velocity gives 



R J V R\ 

Then the null radial acceleration condition ( p2|) gives the well-known result 



12M 2 \ 

Mtf-ftR-ZM): U = R - — ll + yi-— J, (24) 

leading to the requirement R > 6M for stable circular orbits to exist. 

With this in mind, and the explicit formulae for the Christoffel coefficients of 
Schwarzschild metric (given in the Appendix 3), we can establish now the four 
differential equations that must be satisfied by the geodesic deviation 4-vector 
n^- (s) close to a circular orbit. We recall that on the circular orbit of radius R 
(which is a geodesic in the background Schwarzschild metric) we have: 



t dt e dr j, dtp £ d6 

because r = R = const., 9 = n/2 = const., so that sin 9 = 1 and cos 9 = 0. 



4 Geodesic deviation around circular orbit 

It turns out that the four equations are much easier to arrive at if we use the 
explicit form of the first-order deviation equation (|I0|). We get without effort the 
first three equations, for the components n e , and n*: 

V)V = -4n', (26) 



ds 2 y ' R A 
d 2 ^ 2£ dn r d 2 n l 2Me dn 



ds 2 R 3 ds ' ds 2 i? 2 (l - 2M. )2 ds 



(27) 
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The deviation n e is independent of the remaining three variables n*, n r and 



The harmonic oscillator equation (|26| ) for n e displays the frequency which is 



equal to the frequency of the circular motion of the planet itself: 



n e {s) 



n e cos(u s + 7) =n e Q cos (^s + 7^ . 



(2? 



This can be interpreted as the result of a change of the coordinate system, with a 
new z-axis slightly inclined with respect to the original one, so that the plane of 
the orbit does not coincide with the plane z — 0. In this case the deviation from 
the plane will be described by the above solution, i.e. a trigonometric function with 
the period equal to the period of the planetary motion. Being a pure coordinate 
effect, it allows us to eliminate the variable n e by choosing n e = 0. 

It takes a little more time to establish the equation for n r , using Eq. (|i"0|) : 



d 2 n r 



+ 2IV 



0. 



(29) 



ds 2 ' "~ Ap ~ ds 

Taking into account that only the components u f and of the four-velocity on 
the circular orbit are different from zero, and recalling that we have chosen to set 
n e = 0, too, the only non-vanishing terms in the above equation are: 

x dn* 



^ +2 r u t — + 2r 

ds 2 + ttU ds +ZL " 



ds 



+ 9 r r t r t «Vn r + 9J 



0. 



Using the identities 

d 2 rf U 2 



and the definitions 
2M 



5H, we get 



r 2Me dn l 2£ / 2M\ dn^ 
n r + — - — ( 1 - —J 



0. 



ds 2 i? 4 V i2 / _R 2 rfs R V R J ds 

The system of three remaining equations can be expressed in a matrix form: 



(30) 



(31) 



/ ii 

2Me d 



2Me 



V 



R 2 ds 




B 2 (1 _2M)2 ds 

&_ 3i 2 / 1 2M ' 

i? 3 ds 







_2£ /1 _ 2M_\ A. 
R \ R I ds 

ds 2 / 





(32) 



The characteristic equation of the above matrix is 



2M 



4Me 2 



2M\2 
R I 



which after using the identities 



A 2 + 



R ) R\\ 
|) and ( |2^) reduces to 
I 2 ( 6M 



(33) 



1 



R 4 V R 
so that the characteristic circular frequency is 



0. 



CO 



R 2 



1 



6M 
~R 



uj \ 1 - 



6M 
~R' 



(34) 



(35) 



8 



It is obvious that the general solution contains oscillating terms cos(us) ; however, 
before we analyse in detail this part of solution, let us consider the terms linear in 
the variable s or constants: as a matter of fact, because of the presence of first and 
second-order derivatives with respect to s in the matrix operator (p2|), the general 
solution may also contain the following vector: 

(A u*) s + A t \ 

(Au r )s + Ar . (36) 
(Am v )s + Ay?/ 

When inserted into the system (|32|) , the solution is the following: 
A t and A ip are arbitrary; 

Au r = 0, which means that the radial velocity remains null; and 
3£ 2 / 2M\ . 2Mb . . 2£ / 2M\ . w 

This condition coincides with the transformation of the initial circular geodesic of 
radius R to a neighboring one, with radius R + Ar, with the subsequent variations 
A u* and Au^ added to the corresponding components of the 4- velocity in order to 
satisfy the condition u^u v = 1 in the linear approximation. After choosing an 
optimal value for r, we can forget about this particular solution, as well as about 
the arbitrary shift in the variables t and 0, and investigate the oscillating part of 
the solution. 

We shall choose the initial phase to have (with n r > 0): 

n r (s) = — n r cos (us). (38) 

What remains to be done is to compare this frequency with the fundamental cir- 
cular frequency u = £/R 2 of the unperturbed circular orbital motion. 

But this discrepancy between the two circular frequencies u and u is exactly 
what produces the perihelion advance, and its value coincides with the value ob- 
tained in the usual way ([!]) in the limit of quasi-circular orbits, i.e. when e 2 — > 0: 
we get both the correct value and the correct sign. 

Let us display the complete solution for the first-order deviation vector (s) 
which takes into account only the non-trivial degrees of freedom: 

n e = 0, n r (s) = —n r cos(kjs), nf = rift sin(u;s), n* = sin(ws). (39) 

The only independent amplitude is given by rf 0) because we have 



* 2M£ r 2\JM r 

&(l- )2 U ^ 1-2M Jl-Sjf 



- 1 r 2u; r 2 

n o = ~B3~ n o = -B- n o = 1 n o ■ 41 

R A L) RjjJ U _ 6M 



R 
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So the trajectory and the law of motion are given by 



r = R — n T Q cos(c<js), (42) 



(p = ujq s + Uq sin(co>s) = s + Uq sin(a;s), (43) 

R3/2 ./l_3M 



R 



" - /i 2M\ S + ra o sin(ws) = s + n sin(ws), (44) 

where the phase in the argument of the cosine function was chosen so that s = 
corresponds to the perihelion, and s — — to the aphelion. It is important to note 
once again that the coefficient n^, which also fixes the values of the two remaining 
amplitudes, n\ and Hq , defines the size of the actual deviation, so that the ratio ^ 
becomes the dimensionless infinitesimal parameter controlling the approximation 
series with consecutive terms proportional to the consecutives powers of -^f. 

What we see here is the approximation to an elliptic orbital movement as 
described by the presence of an epicycle (exactly like in the Ptolemean system J7], 
except for the fact that the Sun is placed in the center instead of the Earth) . As a 
matter of fact, the development into power series with respect to the eccentricity 
e considered as a small parameter, and truncating all the terms except the linear 
one, leads to the Kepler result [Q, 



r(t) = " (1 6 } ~a[l-e cos(u 1)] , (45) 
1 + e cos(u;o t) 

which looks almost as our formula fl42|) if we identify the eccentricity e with ^ and 
the greater half-axis a with R; but there is also the additional difference, that the 
circular frequency of the epicycle is now slightly lower than the circular frequency 
of the unperturbed circular motion. 

But if the circular frequency is lower, the period is slightly longer: in a linear 
approximation, we have 

I e 2 7 6M\ . . 



R 4 V R , 
hence keeping the terms up to the third order in ^ 



m m ( 3M 27 M 2 135 M 3 \ 

r. ro ^ 1 + _ + __ + __ + ...j. ,47) 

Then obviously one must have ^ = ^r- from which we obtain the perihelion 
advance after one revolution 

6nM IlixM 2 135 ttM 3 . , 

— + + + - (48) 

It is obvious that at this order of approximation we could not keep track of the 
factor (1 — e 2 ) -1 , containing the eccentricity (here replaced by the ratio -£■) only 
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through its square. In contrast, we obtain without effort the coefficients in front 
of terms quadratic or cubic in This shows that our method can be of interest 
when one has to consider the low-eccentricity orbits in the vicinity of very massive 
and compact bodies, having a non-negligible ratio ^. 

In order to include this effect, at least in its approximate form as the factor 
(1 + e 2 ), we must go beyond the first-order deviation equations and investigate the 
solutions of the equations describing the quadratic effects flltt ). 



5 The second-order geodesic deviation 



After inserting the complete solution for the first-order deviation vector (|39|)-(|4l|) 
into the system (^) and a tedious calculation, we find the following set of linear 
equations satisfied by the second-order deviation vector 6 M (s): 



\ 



ds 2 

2Me d 
R 2 ds 





2Me 



H 2 (1 _2M)2 ds 

d? 3^ (i _ 2M_\ 

ds 2 R 4 \ L R ) 
21 d_ 

R' A ds 



2£ 
R 





R > ds 

ds 2 



\ 



J 




n, 



r\2 



(49) 

where we have put into evidence the common factor (^q) 2 , which shows the explicit 
quadratic dependence of the second-order deviation vector 6 M on the first-order 
deviation amplitude n^. The constants C 4 , C r and C v are expressions depending 
on M, R, Uq, lo, e, sm(2us) and cos(2c<js): 



6M 2 (2- 7 -f )esin(2cjs) 



(50) 



3M 



(2- 



5M , 18M 2 ■ 
R R 2 ■ 



-(6- 



27 M 
R 



+ 



6M 2 \ 
R 2 > 



cos (2us) 



2(1 



3M 
R 



)(l-^)i? 4 



6M(1 - f )w sin (2us) 
(1 - 3M )R5 U 



(51) 



(52) 



The solution of the above matrix for ^(s) has the same characteristic equation 
of the matrix (|3^) for n M (s), and the general solution containing oscillating terms 
with angular frequency u is of no interest because it is already accounted for by 
n M (s). But the particular solution includes the terms linear in the proper time s, 
constant ones, and the terms oscillating with angular frequency 2uo: 



Me 



R3{1 _ 6M )(1 _ 2M )2 



3(2 + 



18Af 2 
R 2 



13M 
R 



GM 
R 
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sin(2o;s) 



(53) 



2R(i - m.] 



3(2- 



5M _i_ 18M 2 
R ~r R 2 



1 _ 6M 
1 R 



+ 



5M . 
2 + — ) cos (2us) 
R 



(54) 
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w 



R 2 (l- 



6M\ 



3(2 



5M , 18M 2 ■ 
R R 2 



1 - 



6M 
R 



8M 
R 



2cu 



sin(2u;s) 



(55) 



As explained in Section 2, we need to calculate | 5 2 x' 1 to obtain the geodesic 
curve x' 1 with second-order geodesic deviation: 



5H 



K) 2 ^ 
i? 3 



o/o _ 5M , 18M 2 ^ 
R ' R 2 



S 2 r 



R(l - 



5 2 <p 



6M> 
R > 



_ 33M , 90M 2 
R R 2 



2 - 



15M 
R 



+ 



14M 2 
R 2 



sin(2o;s) 



(i-^)(i-^ 

72M 3 

cos(2a;s) 



R 3 



R 



UJ 



RHl 



6M ' 
R ' 



q/O 5M j 18M 2 
6 \ Z ~ + ~R2- 



6M 
R 



-s H 



R 



32M 
R 



2cu 



sin(2a>s) 



(56) 



(57) 



(58) 



The fact that the second-order deviation vector turns with angular frequency 
2uj enables us to get a better approximation of the elliptic shape of the resulting 
orbit. The trajectory described by x M including second-order deviations is not an 
ellipse, but we can match the perihelion and aphelion distances to see that R^ a 
and e ^ n r Q /R when second-order deviation is used. The perihelion and aphelion 
distances of the Keplerian, i.e., elliptical orbit are a(l — e) and a(l + e). For 
the perihelion is obtained when ujs = 2kir and the aphelion when ujs — (1 + 2k)n, 
where k G Z. Matching the radius for perihelion and aphelion, we obtain the 
semimajor axis a and the eccentricity e of an ellipse that has the same perihelion 
and aphelion distances of the orbit described by 

\2 



R + 



12R 



-1 + 



1 



2M 
R 

6M\2 



+ 



1 



6M 
R 



+ 



15 



(1 



R 



(59) 



2MW-I _ CM 
R Jv 1 R 



n, 



jR (l_2M)(l_6M)2 4 



R 



9M 
R 



+ 



11M 2 , 6M 3 1 
~R2- + "flS-J 



i? 3 



In the limit case of 4r 



0, there is no perihelion advance and a = R 



1 



(60) 

R > 



and e = so the second-order deviation increases the semimajor axis a of a 

Tl r 

matching ellipse compared to the first-order deviation, when a = R and e = 

Another comparison with elliptic orbits concerns the shape of the orbit de- 
scribed by r(ip). From ip(s) it is possible to write s(ip) by means of successive 

approximations, beginning with ujs = ip^J 1 — ^ 

r(s) and we obtain r((p) up to the second order in 



79-1/1 — ^jf. Finally, s can be replaced in 



R 



r 
R 



n' r 



COS 



R 



UJ 



UJq 



-'•r 



3 _ 5M 
R 



30M 2 
R 2 



72M 3 
R 3 



2(1 



2M 
R 



)(1 



R 



2(1-^) 



[2uj 
cos — (p 



+ 



(61) 



(62) 
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In the limit ^ — > 0, the exact equation of a ellipse is obtained up to the second 

[l + e 2 )R: 



order in e, where e = tIq/R and ro 

_ r = (l + e 2 )fl 
1 + e cos ip 1 + e cos ^ 



1 — e cos y> + I - H — cos 2y? 



• (63) 



Comparing with the ellipse equation (|45|), we have r = a(l — e 2 ), so a 
-R(l + 2e 2 ) which agrees with the analysis of Eqs. (^9|)-(l60|). 



R(l+e 2 ) 
(1-e 2 ) 



6 Third-order terms and Poincare's method 

With the third-order approximation we are facing a new problem, arising from the 
presence of resonance terms on the right-hand side. It is easy to see that after 
reducing the expressions on the right-hand side of equation (^) in Appendix 1, 
which contain the terms of the form 

cos 3 us, sin us cos 2 us 

and the like, we shall get not only the terms containing 

sin 3u s , and cos 3 u s 



which do not create any particular problem, but also the resonance terms contain- 
ing the functions sin u s and cos u s , whose circular frequency is the same as the 
eigenvalue of the matrix-operator acting on the left-hand side. 

As a matter of fact, the equation for the covariant third-order deviation can 
be written in matrix form, with principal part linear in the third-order deviation W 1 , 
represented by exactly the same differential operator as in the lower-order deviation 
equations. The right-hand side is separated into two parts, one oscillating with 
frequency u, and another with frequency 3a; : 



d 2 

ds 2 

2Me d 
R 2 ds 





2Me d_ 

jR2(1 _2A L )2 ds 

d 2 3l 2 /-i 2tf ' 

ds 2 R 4 v 1 R ' 

2t d_ 

R 3 ds 







\ 



M (I 

R \ R 



2M \ d_ 
ds 

jL. I 

ds 2 / 




(64) 



r\3 



B f sin(us) + C* sin( 3us) + s D l cos(a;s) 
B r cos{us) + C r cos(3u;s) + s D r sin(a;s) 
sin(cus) + sin( 3us) + s D v cos(us) 



where the coefficients B k , C k and D k , k = t,r,ip are complicated functions of 

The proper frequency of the matrix operator acting on the left-hand side is 
equal to u; the terms containing the triple frequency 3 u will give rise to the unique 
non-singular solution of the same frequency, but the resonance terms of the basic 
frequency on the right-hand side will give rise to secular terms, proportional to 
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s, which is in contradiction with the bounded character of the deviation we have 
supposed from the beginning. The term proportional to s on the right-hand side 
is eliminated in the differential equation for K r when and ^j- are replaced by 
theirs values. 

Poincare pi was first to understand that in order to solve this apparent contra- 
diction, one has to take into account possible perturbation of the basic frequency 
itself, which amounts to the replacement of uj by an infinite series in powers of the 
infinitesimal parameter, which in our case is the eccentricity e = 

uj — > uj + e oj\ + e 2 u>2 + e 3 u 3 + . . . , (65) 

Then, developing both sides into a series of powers of the parameter e, we can not 
only recover the former differential equations for the vectors n M , 6 M , h^, but get 
also some algebraic relations defining the corrections uji, uj 2 , u> 3 , etc. 

The equations resulting from the requirement that all resonant terms on the 
right-hand side be canceled by similar terms on the left-hand side are rather com- 
plicated. We do not attempt to solve them here. However, one easily observes 
that the absence of resonant terms in the second-order deviation equations forces 
Ui to vanish, while the next term uj 2 is different from 0. 

Similarly, as there are no resonant terms in the equations determining the 
fourth-order deviation, because all four-power combinations of sine and cosine 
functions will produce terms oscillating with frequencies 2uj and 4uj; as a result, 
the correction uj 3 will be also equal to 0. Next secular terms will appear at the 
fifth-order approximation, as products of the type cos 5 us, sin 3 ujs cos 2 ujs, etc, pro- 
duce resonant terms again, which will enable us to find the correction uj^, and so 
on, so that the resulting series representing the frequency uj contains only even 
powers of the small parameter 



7 Gravitational radiation 

The decomposition of the elliptic trajectory turning slowly around its focal point 
into a series of epicycles around a circular orbit can also serve for obtaining an 
approximate spectral decomposition of gravitational waves emitted by a celestial 
body moving around a very massive attracting center. 

It is well known that gravitational waves are emitted when the quadrupole mo- 
ment of a mass distribution is different from zero, and the amplitude of the wave 
is proportional to the third derivative of the quadrupole moment with respect to 
time (in the reference system in which the center of mass coincides with the origin 



of the Cartesian basis in three dimensions, see Ref. ||15|| ). 

Of course, it is only a linear approximation, but it takes the main features of the 
gravitational radiation emitted by the system well into account, provided the ve- 
locities and the gravitational fields are not relativistic and the wavelength of grav- 
itational radiation is large compared to the dimensions of the source (quadrupole 
approximation) . 
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More precisely, let us denote the tensor of a given mass distribution \i (xi 



where = 1, 2, 3, see Ref. fTql : 

/ XjdV ^ ^ TTl a X a i X a jj (66) 

where m Q are point masses. 

Let OP be the vector pointing at the observer (placed at the point P), from 
the origin of the coordinate system coinciding with the center of mass of the 
two orbiting bodies whose motion is approximately described by our solution in 
a Fourier series form. It is also supposed that the length of this vector is much 
greater than the characteristic dimensions of the radiating system, i.e. I OP |> R. 

Then the total power of gravitational radiation P emitted by the system over 
all directions is given by the following expression (see Ref. fl6|l ): 

_ G ^d 3 Qij d 3 Qij ^ 1 d?Qa d 3 Qjj\ 

~ 5^ \dl?~ ~dt?~ ~ 3 ~d*3~ ~di?~) ' [ ' 

When applied to Keplerian motion of two masses mi and m^, with orbit equa- 
tion and angular velocity given by 



a(l ~ e 2 ) dip _ y / G(m 1 + m 2 )a(l - e 2 ) 
1 + e cos ip dt r A 

the total power P now reads 

P = 77 — s 2^5 (1 + e cos ^) 4 [12 (1 + e cos (pf + e 2 sin 2 <p] . (69) 
lo c ail. c ) 

We shall calculate the P in Eq. fl67|) with our solution x M using second-order 

M 
R ■ 

. ,, ; / a , can be applied ; 

- 2;? ■ So we finally get P as function of s, which is not shown here because it is 
a very large expression that nevertheless can be easily obtained using a symbolic 
calculus computer program. 

As we want to compare the two total powers P during one orbital period 
(between perihelions) , P in the Kepler case is obtained from the numerical solution 
for (p{t) calculated from Eq. (|68|), and P of the geodesic deviation case has to 
use s{t) obtained from t(s) by means of successive approximations, starting with 

There are many possible ways to compare a Keplerian orbit with a relativistic 
one. Here we assume mi ^> and fix the values of a, e, m\\ the values of R 
and TT-o are calculated to obtain an exact ellipse (up to the second order in e) in 



geodesic deviation, to inspect the non- negligible effects of the ratio We have 
the explicit solutions r(s), tp(s) and t(s), so to calculate ^f- we need only the 
derivatives with respect to s, i.e., % = can ^ e a Pphed successively to obtain 
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Figure 1: The total power P in four cases as function of t during one orbital period T, with 
M = mi. The dotted line is the circular orbit, i.e., e = n 7 = case, and its total power P 
is used as a reference to the others. With e = 0.05 and ^ = 0.04, the total power P for 
elliptic Keplerian orbit is represented by the dot-dashed line. The dashed line is P with x 1 * with 
first-order geodesic deviation, and ^ = 0.05 and ^ = 0.04. Finally, the total power P with 
second-order geodesic deviation is given by the solid line, where ^ = 0.05 and ^ = 0.0402. 

the limit jj- — > 0, like Eq. (|63D , so R = |^~^j a and n r Q = Re. Up to first order 
in e, we have R = a. The choice of M — mi allows the two total powers P to be 
equal when e = and ^ — > 0. Figures 1 and 2 show this comparison for small 
eccentricities and non-negligible ^ ratios. 

Because the emitted total powers P calculated with geodesic deviations depend 
on the ^ ratio, we see that the period is not T = ^=§= (third Kepler's law), but 
an increased one, 



This effect is the direct consequence of the form of angular frequency u that 
appears in the first and higher-order geodesic deviations. 

Another expected feature of Figures 1 and 2 : as e (i.e., ^) is kept small, the 
P using geodesic deviations converge very fast in respect of the orders of geodesic 
deviation. 

Caution is required as the use of quadrupole approximation is not allowed for 
high values of 4r, so the exact amplitude and shape of P using geodesic deviations 




(70) 
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0.2 



Figure 2: The total power P in four cases as function of t during one orbital period T, with 
M = mi. The dotted line is the circular orbit, i.e., e = ?1q = case, and its total power P is used 
as a reference to the others. With e = 0.10 and ^ = 0.02, the total power P for elliptic Keplerian 
orbit is represented by the dot-dashed line. The dashed line is P with with first-order geodesic 
deviation, and ^ = 0.10 and 4| = 0.02. Finally, the total power P with second-order geodesic 
deviation is represented by the solid line, where % = 0.10 and 4! = 0.0201. 



can only be calculated if additional ^ contributions to the gravitational radiation 
formula are included. This approach, but using the post-Newtonian expansion 
scheme, is well developed in Refs. fL7L ELM [19 . 



8 Discussion 

In this paper we have introduced an extension of the geodesic deviation idea in 
order to calculate approximate orbits of point masses in gravitational fields. This 
scheme is of practical applicability to the problem of the emission of gravitational 
radiation. Although in the present paper we restricted our investigations to the 
case of Schwarzschild background fields, our method can be easily extended to 



other background fields [20]. An example is provided by the discussion of Reissner- 



Nordstr0m fields in Ref. J21] 



Since the initial work by Einstein, the problem of orbits and radiation is ad- 



dressed in the literature mostly through the post-Newtonian expansion scheme [16 



2"2|-|]25||. In this approach the starting point for the successive approximations is 
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found in Newtonian theory, whilst relativistic effects are introduced by corrections 
of higher order in - c or The advantage of this approach is, that one can start 
with an orbit of arbitrary high eccentricity. In contrast, our approach starts from 
a true solution of the relativistic problem, but in the case of the Schwarzschild 
background we choose a circular one. We then approach finite eccentricity orbits 
in a fully relativistic scheme, by summing up higher-order geodesic deviations, for 
which we have derived the explicit expressions. 

The two approaches are complementary in the following sense: the post-Newt- 
onian scheme gives better results for small values of ^ and arbitrary eccentricity, 
whereas our scheme is best adapted for small eccentricities, but arbitrary values 
of ^ < |. In both approaches the emission of gravitational radiation is estimated 
using the quadrupole formula, based on a flat-space approximation. 

The next challenge is to include finite-size and radiation back-reaction effects. 
In the post-Newtonian scheme some progress in this direction has already been 
made. In this aspect our result may be regarded as the first term in an expansion 
in ]g. Other applications can be found in problems of gravitational lensing and 
perturbations by gravitational waves. 
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Appendix 1 



The covariant third-order deviation equation is obtained via the same procedure 
that has served to derive the second-order covariant deviation. The third-order 
geodesic deviation itself is 

*V = (p-po) 3 9t 

= e 3 [ft" - artnV + (dsi - 2r»„r;) n " n V] , (7i) 

where h p = Db p / Dp = D 2 n p / Dp 2 . We derive the third-order deviation equation 
by taking the covariant derivative w.r.t. p of Eq. (|T5|) for b p , with the result 



D 2 h p 



Db u 



Ds 2 +R p x> X W= 6R x ^(u x n^- Ds 



+ u' 



~Ds~ 



6") 



n 



"d7 



u v u p n x W + u p u a n v b x 



,Dn a 



.Dn p 



Db° 



+2 I w p n A n^— — + u u n A n a — - + u u u p n x -^- + u p u a n v 



+ 4R.Jn p 



Ds 
Dn x Dn a 



Xp ° Ds Di 



Ds Ds 
+ AR Xp p R a ^u x u^n a n p . 



Db x 
~Ds~ 



(72) 

Related studies of higher-order differentials and their covariant generalizations 



from a more general perspective can be found in recent papers [26, 27 



Appendix 2 

The contributions of various orders to the geodesic deviation obtained in Section 
2 can be deduced in an elegant, coordinate-independent and slightly more general 
manner [2S]. Given a one-parameter congruence of geodesies, one can define the 
tangent vector field Z and the local Jacobi field X; then the Lie bracket of these 
fields vanishes (because the congruence spans a submanifold and therefore is inte- 
grable), so that [X, Z] — 0. 

The geodesic equation is Vz Z = 0. Then, applying the definition of Riemann 
tensor to the vectors X, Y and Z: 

[VxV z - VzV x ] Z - V [x ,z] Z = R(X, Z) Z, (73) 

and taking into account that [X, Z\ = as well as the fact that VxZ = VzX and 
the anti-symmetry of R(X, Z) in its two arguments, we get easily 

V Z V X Z = V 2 Z X = R(Z,X) Z, (74) 
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which coincides with Eq. (|9]) for the geodesic deviation, after we identify the com- 
ponents of the vector fields as Z^ = u M , X M = n/*. 

One may continue in the same spirit and introduce two linearly independent Ja- 
cobi fields, X and Y, both satisfying [X, Z] = = [Y, Z], to obtain the coordinate- 
independent form of Eq. (|T5|), as follows. The two linearly independent Jacobi 
fields, X and Y, satisfy [X, Z] = and [Y, Z\ = 0. By virtue of the Jacobi iden- 
tity, we have [[X, Y], Z] = 0, hence [X, Y] is also a Jacobi field (i.e., it satisfies Eq. 



(|T4|)). Applying the same formula to this field, we get 

V| ([X,Y])=R(Z,[X,Y])Z (75) 

Then, using the fact that Vj Y — Vy X = [X, Y] , we can write the left hand side 
of the above equation as V| (Vx Y — Vy X]). 
Next, we can write this equation explicitly as 

V| (V X Y - VyX) = R(Z, V X Y - VyX) Z 

and furthermore, using the linearity property, as 

V|(VyX) - R(Z, VyX) Z = V 2 Z {V X Y) - R(Z, V X Y) Z. (76) 

The left-hand side of the above equation coincides with the usual Jacobi equation 
applied to the field VyX, whereas the right-hand side can be transformed using 
the definition of the Riemann tensor: the term V|(Vx^) gives 

V|y = Vz (VzVxY) = Vz(VzV x Y - V X V Z Y) + V Z (V X V Z Y) 
= V Z (R(Z,X)Y) + V Z (V X V Z Y) 

= [(V Z R)(Z, X)] Y + R(Z, V Z X)Y + R(Z, X) V Z Y + V Z (V X V Z Y) 

(77) 

(here we used the fact that VzZ = 0). Manipulating further in the same manner 
the commutators of covariant derivations, we arrive at the result 

V| ( Vy X) - R(Z, Vy X) Z = 

Vx R(Z, Y)Z + V Z R(Z, X)Y + 2R(Z, Y)V z X + 2R(Z, X) V z Y (78) 

which is equivalent to Eq. (^) upon identification Y = X and 6 M = (VxX) M . 

Although these coordinate- independent derivations are more elegant, their re- 
sults are not so useful for pratical computations, i.e., the non-manisfestly covariant 
form of the results is better adapted for the calculus of successive deviations in a 
given local coordinate system. 
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Appendix 3 

Connections and curvatures for Schwarzs child geometry 

In this appendix we collect the expressions for the components of the connections 
and Riemann curvature used in the main body of the paper. 

A. Connections. From the line-element ( [L7|) one derives the following expressions 
for the connection coefficients: 



T^a = \ 9*" (d P 9,x + d x g P a - d a9pX )- (79) 



M M ( 2M 

rt L rr _ 2MV U r 2 

r„ = - r(fa ',fl-™Y n. = - r (l-™\ (80) 

r% = K = p ^ = ^. r^ = -sine cose. 

B. Curvature components. The corresponding curvature two-form components 
R P v = 2 Rn\ P udx K A dx^ are: 

Rtr = ^-dtAdr, R te = -f (l - ^) dtAffl, 

M ( 2M\ 2 M 
Rtip = — — [1 — j sin 9dtAd<p, R r e = ^ m ^drAd9, (g;Q 

= ^ 9Mx sin 2 9 dr A dcp, R dip = - (2Mr) sin 2 9d9 A dip. 
r(l — — ) 
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